Numerical convergence of the branching time of negative streamers 
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In sufficiently large gaps and electric fields, discharge streamers do branch. In [Arrayas et a/., 
PRL 88, 174502 (2002)], we observed streamer branching numerically within a deterministic particle 
density model and explained it as a Laplacian instability of a thin space charge layer. Our numerical 
results were criticized in [Kulikovsky, PRL 89, 229401 (2002)]. We here present an adaptive grid 
refinement method for streamer simulations, and we carry out the first conclusive investigation on 
the effect of the numerical grid on streamer branching in different fields. On stepwise finer grids the 
branching time converges, hence streamer branching is for the first time predicted quantitatively. 



Problem setting and review. Streamers are tran- 
sient weakly ionized plasma channels that rapidly grow 
into a non- or weakly ionized medium under influence 
of the self-enhanced electric field at their tip. They 
are widely used in technology [l|, I3 and ubiquitous in 
nature, where the;^play a role in creating the path of 
sparks, lightning 3] and of blue jets above thunder- 
clouds. Streamers are also directly observed as so-called 
sprites 0, IM 5 which are very large discharge structures 
in the higher parts of the atmosphere that are composed 
of ten thousands of streamers. Despite their high ve- 
locity, streamer evolution is now directly observable in 
experiments; a further review can be found in 2]. 

Streamers commonly branch in experiments if gap and 
applied voltage are large enough. Recently a debate has 
risen about the proper physical concept for this branch- 
ing. In 1939, Raether 7] proposed a mechanism for 
streamer propagation and Loeb and Meek |8] developed it 
into a branching concept that nowadays is found in many 
textbooks. The concept is based on a uniformly charged 
streamer head; ahead of it stochastic processes create 
secondary avalanches, that subsequently develop into dif- 
ferent branches. However, the distribution of rare elec- 
trons due to photo-ionization or background ionization 
ahead of the streamer has never been shown to agree with 
the conceptual pictures, and the concept has never been 
demonstrated to work. Furthermore, simulations in the 
past two decades ^, J£., JT, 12] have shown that the fully 
developed streamer head is not homogeneously charged, 
but rather neutral and surrounded by a thin space charge 
layer which enhances the field ahead of it and screens it in 
the interior; this field enhancement allows the streamer 
to penetrate regions with a rather low background field. 
Recent simulations also show that a streamer can branch 
within a fully deterministic model for charged particle 
densities, in a non-uniform background field J3: U3, Um 
as well as in a uniform field [la . UtI Uq , provided certain 
requirements on the external parameters are met (e.g. a 
sufficiently strong background electric field and a suffi- 
ciently long gap). 

Some of the present authors have proposed [l^ [131 
a physical explanation of these numerical observations 
that is directly related to the formation of the thin space 
charge layer: the layer creates an almost equipotential 



streamer head that can undergo a Laplacian instability 
and branch in a manner similar to branching instabil- 
ities of fluid interfaces in viscous fingering. For a fur- 
ther discussion of the conceptual questions of streamer 
branching!:, we refer to B. However, the numerical codes 
used in [ll, [l3, [H [H O El were not able to test the 
branching conditions on fine numerical grids. This lead 
some researchers to question the physical nature of the 
instabilities [l3. llsL [l9l l20J despite the analytical argu- 
ments given in [la:[l3 and later in [211122]. 

To resolve the debate from the numerical side, we 
have developed a code with comoving adaptive grids and 
we here present its results. The algorithm enables us to 
run the simulations on very fine grids; therefore for the 
first time the effect of numerical grids on the branching 
process is investigated quantitatively. We here present 
its results: branching occurs both at very high fields like 
in Refs. [l^ [I^ and also at fairly low background fields 
if the discharge has sufficient space to develop; and the 
branching time saturates on sufficiently fine numerical 
grids. This enables us to give the first quantitative 
predictions on streamer branching. 

Model and multiscale structure of negative 
streamers. We investigate a minimal continuum model 
for streamers, which contains the essential physics for 
negative streamers in a non-attaching pure gas like N2 or 
Ar 9, 10, 16, 17]. The model is a two-fluid approximation 
for the charged particles, with a local field dependent im- 
pact ionization reaction coupled to the Poisson equation 
for electrostatic particle interactions. We investigate this 
model in a cylindrically symmetric geometry, reducing 
it to effectively two dimensions. This constraint sup- 
presses one degree of freedom for the instability modes, 
and therefore the time of branching in this cylindrical 
geometry is an upper bound for the branching time in a 
genuine three dimensional system [2, [23] • In dimension- 
less units, the model reads 
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where a and p are the electron and positive ion densities. 



respectively. E and ^ are respectively the electric field 
and potential, D is the electron diffusion coefficient and r 
is the dimensionless time. The characteristic scales in this 
model depend on the neutral gas density; therefore the 
simulation results can be applied to high altitude sprite 
discharges at low pressures as well as to high pressure 
laboratory experiments. We refer to (2, llS IS^ foi" more 
details on the dimensional analysis. 

A planar cathode is placed at 2: = and a planar an- 
ode at 2: = L^. The potential at the electrodes is fixed, 
0(r, z = 0,r) = 0, (\){r^z = I/2,r) = ^0 > 0, generating 
a background electric field with strength \E\\ = (J)q/Lz 
along the negative z-direction. The streamer is initiated 
by an electrically neutral Gaussian ionization seed on the 
axis of symmetry at the cathode (r = 2; = 0). There is 
no background ionization far from the initial seed. 

We impose homogeneous Neumann conditions for the 
electron density at all boundaries. This results in a net 
inflow of electrons from the cathode if the streamer is 
attached to it jly, |2^. In practice, the computational 
volume is restricted in the radial direction by a bound- 
ary Lr sufficiently far away not to disturb the solution 
near and in the streamer. Moreover, we choose the inter- 
electrode distance Lz so large that the streamer does not 
feel the anode proximity for the results shown. 

The generic spatial structure of the streamer is already 
discussed above and can be seen in the figures: it contains 
a wide range of spatial scales, from the very extended 
non-ionized medium on which the Poisson equation has 
to be solved through the length of the conducting channel 
and its width up to the inner structure of the thin space 
charge layer around the streamer head. 

Moreover, the region just ahead of the streamer where 
the field is substantially enhanced and the electron den- 
sity is low, is highly unstable, in the sense that a small 
ionized perturbation will grow much more rapidly than 
in the mere background field. This unstable region ahead 
of the streamer tip is commonly referred to as leading 
edge 24, 26]. It requires special care when considering 
numerical methods [2^ [2^. Accurate simulations of 
streamers therefore pose a great computational challenge. 

Numerical algorithm. In order to deal efficiently 
with the numerical challenges posed by this model, it 
has been implemented in a numerical code using adap- 
tive grid refinements. We recall the essential features of 
this algorithm and refer to ^] for further details. The 
spatial discretizations are based on finite volumes, using 
a flux limiting scheme to prevent spurious oscillations in 
the results near steep gradients. The time stepping is per- 
formed with a two-stage explicit Runge-Kutta method. 

Using an explicit time-stepping method allows us to 
decouple the computational grids for the continuity equa- 
tions O-iEJ on the one hand from those for the Poisson 
equation (|3|) on the other hand. The particle densities are 
first updated on a series of nested, stepwise refined grids. 
Then the Poisson equation, using the computed densities 
as an input, is solved on another series of nested, step- 
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FIG. 1: Electron density distribution before and just after 
streamer branching in a background field \Eh\ — 0.15, com- 
puted on different finest mesh sizes /i/ = 8, 4 and 2 as in- 
dicated over the plots. The upper snapshots at t= 10000 are 
taken before branching and the lower ones after branching, 
at time r = 11250. The contours correspond to the same 
density levels. In all three cases the same restricted part of 
the total computational domain with z < Lz — 32768 and 
r < Lr — Lz/2 is shown. 



wise refined grids. The electric field on the grids for the 
continuity equations is then calculated from the poten- 
tial computed on the grids for the Poisson equation using 
sufficiently accurate interpolations |27j . 

Adequate refinement criteria for the continuity and 
for the Poisson equation then lead to a grid distribu- 
tion which is especially designed to cope adequately and 
efficiently with the difficulties inherent to both type of 
equations. More specifically, the refinement criterion for 
the grids for the Poisson equation is based on error esti- 
mate of the solution. The refinement criterion for grids 
for the continuity equations uses a curvature monitor of 
the solution. Moreover, it takes explicitly into account 
the leading edge^ where the densities are low but the elec- 
tric field is greatly enhanced [2J, [2y| . 

The refinement criterion is computed at each time 
step, in such a way that the series of nested, consec- 
utively refined grids move with the solution. Special 
care has been taken for the discretizations as well as the 
mapping of the solution from one grid to the other to be 
charge conserving. 

Results. The adaptive grid refinement procedure en- 
ables us to resolve the streamer with very high accuracy, 
and thus to investigate the dependence of the branching 




T=10250 



h=1/2 



600 



400 



200 



600 



h=1/4 h^=1/8 

600 1— ■ — i600r 



400 400 400 

# .00 # .00 ® .00 ® 

-100 100 -100 100 -100 100 



I 



0.6 
0.4 
0.2 




600 



400 



200 




100 100 

r 



T=10750 T=11250 

900 F 




100 200 100 200 

r r 

T=325 T=350 

550g =~^^^^v^v.v. i 600 1 



100 200 



r 

T=375 




50 100 150 50 100 150 



50 100 150 



FIG. 3: Zoom into the streamer head during branching. Up- 
per plots: \Sb\ =0.15 as in Fig. Q /i/ = 2. Lower plots: 
\Sb\ =0.5 as in Fig. |21 hf = 1/8. Contour lines (thick) of net 
charge density and equipotential lines (thin) are shown as a 
function of positive radius r and appropriate z. The spacing 
of the charge contour levels is 0.004 for the low field case, and 
0.16 in the high field case. The spacing of equipotential lines 
is 5 in both cases. 



FIG. 2: Upper panel: Branching time in a background field 
\Sb\ = 0.5 as function of the finest mesh size hf = 2, 1, 1/2, 
1/4, 1/8. Lower panels: the corresponding electron density 
distribution at r=275 (middle row), and just after the respec- 
tive branching time (lower row), computed on different finest 
grids hf = 1, 1/2, 1/4 and 1/8. The total computational 
domain is 2; < L^ = 2048 and r < Lr — Lz/2. 



process on the numerical grid. The results are obtained 
on increasingly finer grid sizes /i/, always taking the 
same coarsest mesh width he for both the continuity 
and the Poisson equations. If the branching were of 
numerical nature, we would expect that branching times 
on increasingly finer grids would not converge. 

We first consider negative streamers evolving in a low 
background field of S^ = 0.15 corresponding to 30 kV/cm 
for N2 at atmospheric pressure. We use an electrically 
neutral, dense and relatively wide Gaussian ionization 
seed at the cathode, with a maximum of 1/4.8, and a 
characteristic radius of 10. This corresponds to a max- 
imal electron and ion number density of 10^^ cm~^ and 
an 1/e radius of 230 /im. The gap length and width are 
set to Lz = 2Lr = 2^^ = 32768, which corresponds to an 
inter-electrode distance of approximately 7.5 cm. 

The coarsest mesh width is set to he = 64, and the 
finest one to /i/ = 8, 4 and 2. When a finest mesh of 8 
is used, the electron density in the streamer is lower, as 
can be seen in the upper row in Fig. ^ This is due to 



the numerical diffusion introduced on such a coarse grid 
by the flux limiter that switches to the diffusive first or- 
der scheme in regions with large gradients. This numer- 
ical diffusion smears the electrons out over the streamer 
head, which in turn results in lower field enhancement 
and lower ionization rates. The results on finer meshes 
of 4 and 2 on the other hand do agree with each other. 

The branching in time is the same in all cases. Fig. Q 
shows that the influence of the numerical grid on the 
branching state rapidly decreases, and we thus can carry 
out not only qualitative but also quantitative numerical 
experiments of the streamer evolution up to branching. 
These results show that branc hing is possible at lower 
electric fields than those of [T^ [13 • Branching was not 
observed in earlier simulations at lower fields |3, llfl 
because the discharge gap was too short. 

We now consider a negative streamer in a dimen- 
sionless background field of S^ = 0.5 corresponding to 
100 kV/cm in N2 at atmospheric pressure in a gap of 
Lz = 2048, or 4.6 mm. These external parameters are 
as in 16., J/7j . The initial seed is also taken as in U^ , 
i.e., a Gaussian with amplitude 10~^ and characteristic 
radius 10, which corresponds to a maximal electron and 
ion number density of 5 • 10^^ cm~^ and an 1/e-radius of 
23 /im for N2 under normal conditions. 

However, while 16] used a uniform grid of h = 2 and 
[13 one of h = 1, we now perform computations on a 
finest grid as small as hf = 1/8, i.e., more than a decade 
finer. More precisely, the coarsest mesh width is set to 



he = 2, and the finest one to hf = 2, 1, . . . , 1/8. Further- 
more a better numerical scheme is used: flux limiting [23 
rather than 3rd order upwind 16, 17]. 

Before branching, at r=275, Fig.Oshows that there is 
a quantitative difference between the results on a mesh 
with hf = l and the other three. As in the low field case, 
numerical diffusion spreads the space charge layer, which 
makes the field enhancement at the streamer tip and 
the field screening in the streamer body less efficient. 
Consequently, the ionization rate, and therefore the 
electron density, are higher in the streamer body. In 
the low field case we do not observe this because the 
background field is negligibly low, hence a less efficient 
screening will not affect significantly the ionization rate 
in the streamer body. It is clear that on meshes finer 
than 1/2, the results are the same during the stable 
streamer propagation. It is only after the branching that 
different states are observed on those very fine grids. 
However, the time of branching converges within this 
range of mesh widths hf as shown in the upper plot in 
Fig. 121 

Discussion, conclusion and outlook. We empha- 
size that the branching times converge on decreasing nu- 
merical grids in both cases. Therefore we here present 
the first conclusive and quantitative numerical predic- 
tions on streamer branching. However, in contrast to the 
low field case, the lower plots in Fig. [2 show that in the 
high field case different branched modes are reached af- 
ter approximately the same evolution time: in two cases. 



the maximal electron density and field is on the axis of 
symmetry, and in two other cases, it is off axis. Ap- 
parently, there are different branched states reachable at 
bifurcation and tiny differences determine which one will 
be reached. Such extreme sensitivity is well-known from 
deterministic chaos; it is generic for nonlinear dynamics 
near bifurcation points. On the other hand, the unstable 
state is reached in a deterministic manner, and therefore 
the branching times converge. 

But why is there once a unique branched state and 
once several? The answer can be found in Fig. |31 showing 
the two relevant spatial scales, namely the thickness of 
the space charge layer and the radius of the channel. In 
the high field case, the ratio of layer thickness over radius 
is much smaller than in the low field case. Moreover, the 
field screening and enhancement is much stronger and 
the equipotential lines follow the space charge layer much 
better. Therefore the high field streamer is much closer 
to interfacial models as discussed in y, Ua? IH 123, 124^ and 
can access more branching modes. This critical state in 
future work will be characterized by the electric charge 
content and electric field and potential at the streamer 
tip which would then allow us to relate branching to the 
external electric circuit. For sketches of such ideas as 
well as for a discussion of photo-ionization effects and of 
continuum versus particle models, we refer to 2]. 
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